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Abstract 



Context. New observational means such as the space missions CoRoT and Kepler and ground-based networks are and will 
be collecting stellar pulsation data with unprecedented accuracy. A significant fraction of the stars in which pulsations 
are observed are rotating rapidly. 

Aims. Our aim is to characterise pulsation modes in rapidly rotating stellar models so as to be able to interpret 
asteroseismic data from such stars. 

Methods. The pulsation code developed in Lignieres et al. (2006) and Reese et al. (2006) is applied to stellar models 
based on the self-consistent field (SCF) method (Jackson et al. 2004, 2005, MacGregor et al. 2007). 
Results. Pulsation modes in SCF models follow a similar behaviour to those in uniformly rotating polytropic models, 
provided that the rotation profile is not too differential. Pulsation modes fall into different categories, the three main 
ones being island, chaotic, and whispering gallery modes, which are rotating counterparts to modes with low, medium, 
and high Z — |m| values, respectively. The frequencies of the island modes follow an asymptotic pattern quite similar 
to what was found for polytropic models. Extending this asymptotic formula to higher azimuthal orders reveals more 
subtle behaviour as a function of m and provides a first estimate of the average advection of pulsation modes by 
rotation. Further calculations based on a variational principle confirm this estimate and provide rotation kernels that 
could be used in inversion methods. When the rotation profile becomes highly differential, it becomes more and more 
difficult to find island and whispering gallery modes at low azimuthal orders. At high azimuthal orders, whispering 
gallery modes, and in some cases island modes, reappear. 

Conclusions. The asymptotic formula found for frequencies of island modes can potentially serve as the basis of a mode 
identification scheme in rapidly rotating stars when the rotation profile is not too differential. 



1. Introduction 



New observational means are and will be collecting stellar pulsation data with unprecedented accuracy. The space mission 
CoRoT has considerably lowered the detection threshold for pulsation modes, thus allowing photometric observation of 
solar-like pulsations in stars other than the Sun and increasing the number of detected modes in early- type stars. The 
forthcoming space mission Kepler will add a wealth of pulsation data by observing a large number of stars for a period 
of four years. Other projects include the space mission PLATO as well as ground-based networks such as SONG. 

Stellar pulsations yield valuable information on the internal structure of stars which can be used to constrain stellar 
evolution models. Although a great deal of success has been achieved in probing the internal structure of the Sun 
and of a number of other stars, a number of difficulties arise for rapidly rotating stars. Indeed, rapid stellar rotation 
introduces a number of phenomena which considerably complicate their modelling and the study of their pulsation modes. 
These include cen trifugal deformat ion, gravity darkening, baroclinic flows and various forms of turbulence and transport 
phenomena {e.g. |Rieutord|| 2006b[ ) . As a result, the internal structure of these stars remain difficult to probe. 

Traditionally, the effects of rotation on pulsation modes have been modelled using the perturbative approach. In 
this approach, rotation is taken into account through corrections which are added to the non-rotating solutions. The 
underlying assumption in this method is that the rotation rate, f2, can be treated as a small parameter, th us enabling 
one to develop the perturb ative cor rections as a power series in Such series can be extended to fir st ( Cowling & 
Ledoux|[l95l]), second (|Sa io||1981j |Gough fc Tho mpson fl990[ |Dziembowski fc Goode||1992[ ) or third order 
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oufi et al. 1998 J^v^arami et a!l.' 2005). A natural question to ask is up to what rotation rate is this approach 



valid, i'his remained an open question until ,Reese et al. ( 2006, ) applied a non-perturbative two-dimensional approach 



to calculating acoustic pulsations in polytropic stellar models and compared the results with perturbative calculations. 
Their results showed that perturbative methods remain valid only for values which are l ower than the rotation rat e of 
many early-type stars. Further comparisons between the two appr oaches include those by Lovekin & Deupree (20081, in 
which more realistic models were used but at a lower accuracy, and Ouazzani et al. (2009 1, in which the effects of avoided 
crossings are included in the perturbative calculations. 
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Due to the limitations of the perturbative method, a number of recent studies have focused on modelling the effects 



of rapid rotation on stellar acoustic pulsations using a two-dimensional approach. Espinosa et al. (2004) studied the 



effects of rapid rotation on frequency multiplets in models with a uniform density and also briefly discussed pulsations of 
realistic models. They showed how rotation leads to highly non-uniform multiplets and causes th e frequencies of adjacent 



modes to pair up, thus providing a tentative explanation for observed close frequency pairs (Breger & Pamyatnykh 



20061. Lignieres et al. (20011 studied pulsation modes in a uniform density spheroid using a perturbative method and 



two different numerical approaches. This was done in order to va lidate their two-dimen s ional numerical me thod before 
apply i ng it t o more realistic models. Their work was followed by Lignieres et al. (20061, Reese et al. (20061 and Reese] 



et al. (2008a I who did the first accurate calculations of p-modes in rapidly rotating polytropic models. They investigated 
the limits of the perturbative approach, studied disk averaging factors which intervene in mode visibility, compared the 
effects of the centrifugal and Coriolis forces and found an empirical formula which charac terises the structure of the 
frequency spectrum for low degree modes. At the same time, Lignieres & Georgeot (2008) and Lignieres & Georgeot 
(2009) applied ray dynamics to the study of acoustic modes in rotating polytropic models. They classified modes into 



several categories, the main ones being island, chaotic, and whispering gallery modes which are rotating counterparts 
to modes with low, intermediate, and high I — \m\ values, where I is the harmonic degree and m the azimuthal order. 
They showed that each category has its ow n frequency organisat ion a nd provided an expl anation involving t ravel time 
integrals for the empirical formula found in 
( 12008 J and ,Lovekin et al. 
Deu 



Lignieres et ah (12006') andlReese et al.l (|2008a|). Finally Lovekin 



^ D eupree 

( p009j ) studied p-modes with low radial orders in realistic models from ^Deupree, ( ^1990 ) and 



eupree ( 1995 ) with both uniform and differential rotation. They investigated how frequencies and the large and small 



separations vary with uniform or differential rotation and compared their calculations with a perturbative approach. 

Before being able to interpret pulsation modes in observed stars, more progress is needed in understanding the effects 
of rotation on pulsation modes. Indeed, although a number of important results have been established for p-modes in 
polytropic models, these need to be extended to more realistic models. The calculations involving more realistic models 
have currently been limited to small mode sets and the analysis has not been pushed far enough to see whether similar 
results apply. In what follows, we calculate pulsation modes, using the numerical method developed in 'Lignieres et alT] 
( 2006 ) a nd Reese et al. J 2006 ), in realistic models of rapidly rotating stars based on the Self-Consistent Field (SCF) 
method ( | Jackson et al.||2005 MacGregor et al.|2007 ) . In particular, we investigate whether a similar mode classification 
exists in these models, whether a similar empirical formula applies to frequencies of modes with low ^ — |to| values, and 
quantify the effects of using a differential profile. The next section deals with the SCF method and the models it produces. 
The following section explains the pulsation equations, the numerical method used for calculating the pulsation modes 
and a number of tests to validate the method. Afterwards, Sections |3] and |4] describe the results for models with mildly 
and strongly differential rotation, respectively. Our conclusions are summarised in Section [5j 



2. Stellar models based on the SCF method 

The SCF method is an iterative procedure for solving the equations that govern the structure of a conservatively rotating 
star. The basic approach underlying the method is to alternately solve Poisson's equation to derive the 2D shapes of 
equipotential surfaces, and the equations of mass, momentum, and energy conservation to obtain the ID profiles of 
thermodynamic quantities along a radius in the rotational equatorial plane. As described in detail in Jackson et al. 
(2005), this procedure yields a sequence of models which, under most circumstances, converges to a model that satisfies 
all the equations for a prescribed internal rotation law. 

The metho d was first developed a nd used 40 years ago to compute uniformly and differentially rotating polytropic 
stellar m odels ( jOstriker fc Mark|1968 ) . Although subsequently extended through the incorporation of more realistic input 
physics ( jJackson|1976[ ), application of the method was limited to massive stars, a consequence of conver gence difficultie s 
encountered in lower niass models with sufficiently high values of the central mass concentration (see, e.g., Clement|1978 ). 
This problem was addressed and remedied through a reformulation of the method in which the normalised distributions 
of thermodynamic quantities and the central values of those quantities are adjusted in separate iterative loops. The new 
SCF method has been implemented in a code th at utilises up-to -date input physics. The opacities are obtained from the 
ta bles of OPAL opacities compu ted by Rogers & Iglesias ( 1992 ) and from ta bles of low-tempera ture opacities compiled 
by Alexander & Ferguson (1994), using interpolation subrouti nes written by Vandenberg (1983). The equation of state 
for the stellar material is calculated according to the formula of Eggleton et al. ( |1973 ), and the nuclear energy generation 



Graboske et al. 



rates for hydrog en bu rning are from Caughlin & Fo wler] ( 1988] ), with the treatment of electron screening effects from 
1973| for the case of equilibrium abundances of CNO isotopes. Energy transport in sub -photospheric 



convective envelopes is treated using a standard mixing-length model (see, e.g., Kippenhahn et al. 1967), in which the 
local gravitational acceleration g is replaced by the value of g as reduced by the local centrifugal acceleration, averaged 
over equipotential surfaces. For the models utilised in the pulsation mode computations described in subsequent sections, 
a value of 1.9 was adopted for the ratio of the mixing length to the pressure scale height. The method and the code 
are both robust and rapidly convergent, and have been thoroughly tested and validated thro ugh such application s as 
the interpretation of interferometric observations of rapid rotators like the Be star Achernar ( jJackson et al. 2004 and 
refer ences therein) and an ex amination of the effects of differential rotation on the structure of stars less massive than 2 
Mq |MacGregor erar]]2007t . 

MoHels computed using the SCF method are chemically homogeneous, ZAMS models with the following abundance 
fractions by weight of H, He, and heavy elements: X = 0.7112, Y = 0.27, and Z — 0.0188. The rotation proffie 
is imposed beforehand and is conservative, i.e. the centrifugal force derives from a potential. As a result, the stellar 
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structure is barotropic - different thermodynamic quantities remain constant along lines of constant total (centrifugal 
plus gravitational) potential. The rotation profile used in the present calculations was: 



as 



(1) 



where s is the distance from the rotation axis, R^q the equatorial radius, and flcr the break-up rotation rate at R^q. The 
parameters rj and a determine how rapid and differential the rotation is. fn particular, the ratio between the polar and 



equatorial rotation rate is 1 



The associated angular momentum increases with s thereby satisfying the dynamical 



part of the Solberg-H0iland criterion for stability. Various forms of she ar instabili ty may, however, be p resent if the 
rot ation profile becom es too differential, i.e. if a becomes too large {e.g. Zahn|[l974 ) . Also, as explained in Zahn (19931 



and Rieutord| ( [200 6a') , baroclinic flows occur in radiative zones of rapidly rotating stars thus leading to a non-conservative 
rotation profile. Exploring these effects is, however, beyond the scope of this paper. 

Equation ([!]) corresponds to a rotation profile in which the rotation rate decreases with s. Such profiles can be used 
to construct highly distorted configurations. Indeed, the stellar core can be made to rotate quite rapidly since the local 
break-up v elocity is larger than at the equator. This type of model was used to try to explain Achernar's extreme 



oblateness (Jackson et al. 20041. The SCF method can also produce models with a rotation rate that increases with 



distance from the rotation axis. This resembles somewhat the solar rotation profile in wh ich the rotation rate increases 
with decreasing latitude in the convection zone ( Schou et al.||1998 Thompson et al.||2003 l. 



2.1. The pulsation equations 

In order to derive the set of equations which govern acoustic pulsation modes in a differentially rotating star, we start by 
representing the differential rotation by a permanent background flow Vq = sfl{s)e^. In what follows, we will work with 
cylindrical coordinates {s,z,(j)) and their associated unit vectors (bs, Sz, e^). We write out the Eulerian perturbation to 
various equations, starting with Euler's equation, and only keep first order linear terms: 



Po 



dv 

dt 



pVo ■ VVo + PoV ■ VDo + PoVo ■ Vt> = - Vp + pQo - Po^^ , 



(2) 



where quantities with the subscript "o" are equilibrium quantities, and those without a subscript Eulerian perturbations. 
The quantity go is the background gravity excluding the centrifugal acceleration. The different terms on the left hand 
side of Eq. ([2| can be worked out explicitly in terms of ri(s): 



Vo ■ Vt> — n X V + imUv, 



(3) 
(4) 
(5) 



where we have assumed an e*""^ azimuthal dependence for v and — ile^. The first term corresponds to the centrifugal 
acceleration and the sum of the next two includes the Coriolis force. Combining these equations with Eq. ^ yields: 



(6) 



Oil 

[X + imH.] poV = -Vp + pQctf - PoV^' - 2$! X poV - poS — Vge^ 

where we have assumed an e^* time dependence for v and where gcs — ^Po/po is the background effective gravity which 
includes the centrifugal acceleration. The Eulerian perturbation to the continuity equation gives: 

dp 
dt 

In terms of fi, this becomes: 



V • (pov) + V • (pVo) 



0. 



(7) 



[A -I- imfl] p —V ■ Vpo — Po^ ■ V. 

The Eulerian perturbation to Poisson's equation is simply: 



(8) 



(9) 



where G is the gravitational constant. These equations are then supplemented by the adiabatic relation between the 
pressure and density perturbations which takes on the following form: 



dp 
dt 



V ■ Vpo 



Vp - cl 



dp 
dt 



V ■ Vpo 



Vn 



Vp 



(10) 



where = ^Jj^ is the square of the sound velocity and Fi the adiabatic exponent. Provided that Vo ■ Vpo = Vg ■ Vpo = 
Vo ■ = 0, this form can be shown to be equivalent to Sp — c^Sp where 6p and Sp are the Lagrangian pressure and 
density perturbations, respectively. This leads to the following equation after some manipulations: 



[A + imn] {p - cIp) = [-Vpo + clVpo] ■ v 



(11) 
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2.2. Non-dimensional form 

These equations are then put into non-dimensional form using the following length, density and pressure scale factors: 
Rcq, Pc, Pc, (12) 



where the subscript "c" refers to the centre of the star, and i?eq is the equatorial radius. This gives a time scale t 




'ref 



(13) 



Based on these scale factors, all of the above equations remain the same as in dimensional form, except for Poisson's 
equation where a non-dimensional factor A — AnGp^R^^/pc appears: 



A^f = Ap. 



(14) 



2.3. Splieroidal geometry 

In order to achieve higher accuracy when solving these equations numerically, a coordinate system which follows the 
shape of the star is introduced. This new coordinate system ((^, 9, (f) can be related to the usual spherical coordinate 
system (r, 6, 4>) via the following relationship: 



r(C, 0) = (1 - e)C + ^-^-^ {R,{e) -l + e), 



(15) 



for C g [0, 1]. When C, = 1, r coincides with the stellar surface, Rs{6). The variables 9 and remain the same in both 
systems. A second domain is added around the first, in which r is given by: 



r(C, e)^2e+{l- £)C + (2C' - 9C' + 12C - 4) {R,{e) - I - e) . 



(16) 



for C e [1,2]. With these definitions, r and rc^ = d(^r remain continuous across C, — 1. As C, approaches or 2, this 
coordinate system behaves like a spherical coordinate system: the constant C-lines become spherical and becomes 
independent of 6. This is important as it simplifies the regularity conditions in the centre and the boundary condition 
on th e pertu rbat ion to the gravity potential on the outer boundary. The sa me coordinate system was used in |Ligniere"i| 
et al. (2006) and Reese et al. (20061 and is based on Bonazzola et al. (1998). The stellar model is then interpolated onto 
a grid based on this new coordinate system. 

Another possibility would be to base the radial coordinate on the equipotentials. This has the advantage of simplifying 
the pulsation equations because terms such as dgpo and dgPo vanish. However, this requires using numerical rather than 
analytical differentiation when calculating terms with radial derivatives such as tq , thereby reducing the accuracy of the 
results. Furthermore, the regularity conditions in the centre of the star become more complicated as the equipotentials 
do not in general become circular towards the centre. 

Based on the coordinate system presented above, the continuity equation becomes: 



rx , ■ oi C^cPo c ^^oPo e C Po 
[X + imll\p^ — u'' z^^u - 



dg (sin 9u^ 
Csin6' 



+ 



Csin0 



(17) 



Euler's equation takes on the following form: 



[X + imfl] Po 



[X + imfi] Po 



the adiabatic relation is given by 

e 



Po- 



2nC, sin Bu'^ 



-a, 



xp- 



'jQPo 
Po 



Po- 



2nCirg sin 61 + r COS 6*) m"^ 
2nC^ sin Ou'^ 



p - Pod(;^, 



dgP' 



Po 



-P- PoOg^, 



-Po- 



Po- 



2n({rgsm9 + rcos9)u^ d^p 



-Po-T 



sin 9 



[X + im^] [p - cIp) 



Po sin 9 (dgfl) 



c 



C^sm9u'^ 



C{rg sin ( 



-t- r cos ( 



{~9^Po + cl^CPo) u'^ + {-dgPo + cldgpo) u' 



(18) 
(19) 

(20) 
(21) 
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and Poisson's equation takes on the following form: 



r + r Ira 1 



(22) 



where m^, u^, and v!^ are the th ree veloci ty c omponents (see below for details) and terms of the form r^, tq ... are different 
derivatives of r based on Eqs. (15 1 and (16 1. Explicit expression for and A^^ are as follows: 



[irc^rgrc^B - t^tqc, - rlreg + 2rr^ - rjr^t: ~ r^ro cot ( 



9L + cot Qds 



-dl 



(23) 
(24) 



The above expressions are obtained by using tensorial expressions for the differential operators which intervene and 
working them out explicitly. Furthermore, the components to the velocity are written on the following basis: 



2 



Ea 



{rgSr + reg) 



(25) 



r^r^ sin ( 



where {E^} and {6^} are the natural and spherical basis, respectively. When the star becomes spherical, the basis {ai} 
converges to jeA Apart from some multiplicative factors, Euler's equation is expressed on the dual basis, as is done in 
Reese et al. ( 2006[ ), since this has the advantage of limiting the effective grav ity to the radial component at the stellar 



surface. More details on how to derive the above expressions can be found in Reese ( 2006 1 and references therein 



Equations (17l-(22l were completed with a number of boundary conditions which ensure that the solution remains 



regular in the centre, the Lagrangian pressure perturbation vanishes on the stellar surface and the perturbation to the 
gravitational potential goes to zero at an infinite distance from the star. 



2.4. Numerics 



The above equation s were th en projected onto the spherical harmonic basis in the same way as was done in Lignieres 
et al. ( 2006 ) and fieese et al. ( 2006) . This is achieved by expressing the different unknowns as a sum of scalar or vectorial 
spherical harmonics multiplied by unknown radial functions, and then projecting the equations themselves onto the 
spherical harmonic basis, using Gaussian quadrature to numerically perform the integrations. The resultant system is an 
infinite system of coupled ordinary differential equations in terms of the radial variable C which is truncated at a maximal 
harmonic degree Lmax- The solution to this system yields the radial functions used in the harmonic decomposition of the 
different unknowns. 

This system of ordinary differential equations is discretised using one of three methods: a spectral method based on 
Chebyshev polynomials, finite differences or a polynomial spline-based method. In the latter two cases, the order can be 
adjusted. When applying these methods, the stellar model is interpolated onto either a higher resolution uniform grid or 
a Chebyshev Gauss- Lobatto collocation grid usi ng cubic spline interpolation. For the spectral method, this is analogous 
to what was done in Dintrans & Rieutord (2000 1 in which a 1.5 Mq CESAM model was interpolated onto the same type 
of collocation grid before being used to calculate gravito-inertial modes. 

After discretisation, the eigenvalue problem is in algebraic form: Av = XBv where A and B are square matrices. With 
a suitable choice of variables, these matrices can be made real, thereby reducing the computational cost. Also, pulsation 
modes are either symmetric or anti-symmetric with respect to the equator, so that only spherical harmonics of the same 
parity are need ed to describe them. The problem is then solved numerically using the Arnoldi-Chebyshev algorithm {e.g. 
Chatelin|1988 1 around different target frequencies, called frequency shifts. In what follows, most of the calculations have 
been done using a 4**^ order finite difference approach. The angular resolution was typically L^^x = 80 and the radial 
resolution = 501. 



2.5. Accuracy of the calculations 

Various tests can be used to assess the accuracy of the calculations. A first test consists in following the evolution of the 
frequency error as a function of the radial and angular resolution. The solid lines in Fig. [T] show the evolution of the 
relative error on the numerical frequency for two modes in a 25 Mq model uniformly rotating at 60% of the break-up 
rotation rate, using the frequencies calculated at highest resolution as references. The first two panels apply to an ri = 16 
mode and the other two to n = 50 (see Section for the definition of h). As is evident from the figure, the stability 
of the numerical frequencies is very good, especially for the angular resolution where spectral convergence seems to be 
achieved. 
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Intermediate radial order High radial order 




Figure 1. Evolution of various forms of the frequency error with L^ax and N,- for two pulsation modes in an 25 Mq star 
rotating uniformly at 60% of the break-up velocity. The mode on the left corresponds to n = 16 and the one the right 
to n = 50 (the meaning of n is given in Section |3] and illustrated in Fig. [s]) . The solid and dotted lines correspond to the 
relative frequency error. This first case uses the numerical frequency and the second is based on the variational frequency. 
The dashed curve shows the relative difference between the numerical and variational frequencies. The frequencies fi^ 
and ri™'' are the numerical and variational frequencies calculated at highest resolution {i.e. Lmax = 100, — 751 for 
the mode on the left and Lmax — 160, A^r — 1001 for the mode on the right). As can be seen in the figures, the numerical 
frequencies are very stable as a function of the resolutions, and the variational frequency somewhat less stable. Also, a 
discrepancy remains between the two types of frequencies. 



The right two panels of Fig. 2 of iReese et al. (2008b I show similar curves for a pulsation mode in a 1.8 star 



rotati n g unifo rmly at 90% of the break-up rotation rate. In this case the results were not as good. As explained in Reese 



et al. (2008b I, evaluating the error in this case was not entirely straightforward due to difficulties in identifying the 
correct mode at different resolutions. Indeed, at such high rotation rates, regular modes interact much more with chaotic 
ones thus distorting their geometric features. Furthermore, the amount of interaction between the different modes seems 
to depend on the numerical resolution. 

Although the problem is expressed in terms of real matrices and frequencies are searched for around real target 
frequencies, complex conjugate solutions sometimes appear. For instance, two of the calculations in the panel to the far 
right of Fig. [T| (at = 601 and = 801) correspond to complex solutions. The imaginary parts are most likely due 
to numerical inaccuracies as these solutions are replaced by real solutions at other numerical resolutions. Their relative 
magnitude (3 x 10~^ — 10~^) suggests a comparable accuracy on the corresponding frequencies. 

Another test consists in applying a variational formula on the eigenmodes to yield an independent value for the 
frequency. According to the variational principle, the error on the "variational frequency" is proportional to the square of 



the error on the eigenmode, thus minimising its effect provided it is sufficiently small (Christensen-Dalsgaard & Mullan 



19941. By comparing this value to the original (numerical) frequency, it is possible to estimate the accuracy of the 



calculation. In what follows, we calculated variational frequencies using the following formula which is only valid for 
uniform rotation: 

= {uj^e,r + mnf / po\\vfdV + 2i{uj^^, + mn) / ■ {v* x v) dV - / poN^ \v ■ eg\^ dV 
Jv Jv Jv 



\Jv PoCi J^Jv^ 



(26) 



where uj is the numerical frequency, tJvar the variational frequency, m the azimuthal order, V the volume of the star, Voo 
infinite space, and Sg the unit vector in the same direction as the effective gravity. The geometric term mfi comes from 
the fact that the pulsation frequencies are expressed in an inertial frame. In Section [3. 4| we give a more general variational 
formula which is also valid for differential rotation, but is expressed in terms of the Lagrangian displacement rather than 
the Eulerian velocity perturbation. Such a formulation gives comparable results as Eq. 26 i.e. 5u}/u} < 10^"^ — 10^^, even 
for the most differential rotation profiles. 

The dashed lines in Fig. [ijshow the relative difference between the numerical and variational frequencies. As can be 
seen in the figure, these differences are much larger than the variations caused by modifying the resolution. A third set 
of curves, the dotted lines, show the relative error on the variational frequencies when using the variational frequency at 
highest resolution as a reference. From these, we deduce that the variational frequencies do converge to a specific value, 
but at a slower rate than the numerical frequencies, which is the opposite of what we would expect from the variational 
principle. Furthermore, the limit of the variational frequencies is different than that of the numerical frequencies, as can 
be seen from the dashed curves. Such a discrepancy can occur if the models are not in perfect hydrostatic equilibrium due 
typically to numerical inaccuracies. Indeed, hydrostatic equilibrium is implicitly assumed when deriving the variational 
formula, and deviations from this state will tend to produce errors which are independent of the resolution of the 
eigenfunctions. Nonetheless, these discrepancies remain small when the rotation rate is not too close to break-up and 
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probably affect tfie different modes in a similar way for a given model so that the analysis in the rest of the article is not 
likely to be affected. 




M = a5.0Vt, n o.e a 0.0 
U-279.Bji.m m = 0" 

Finite differences 



U = ^b.OSU fi 0.8 a 0.0 
cj-279.6/iIlz m = 0" 

Chebyshev polynomials 



M=25.0M,, n O.S a 
U- 279.5/.IR7. m = 0" 

Polynomial splines 



Figure 2. Comparison of a pulsation mode calculated using finite differences (left), a spectral method (centre), and 
polynomial splines (right) to discretise the equations in the radial direction. These and other similar plots show a 
meridional cross-section of the Eulerian pressure perturbation divided by the square root of the background density 
profile so as to bring out near surface regions. The difference on the frequencies is less than 0.1 /iHz. 



Finally, a last test consists in applying different numerical techniques to calculate the eigenmodes and seeing if they 
give similar results. Figure [2] shows such a comparison. The mode on the left is calculated using 4'^ order finite differences 
in the radial direction, the one in the middle a spectral method based on Chebyshev polynomials and the one on the 
right 4**^ order polynomial splines. The angular resolution was very similar for the three cases and the radial resolution 
went from = fOl for the calculation based on Chebyshev polynomials to = 30f for the spline-based calculation. 
As can be seen in the figure, the three calculations yield very similar results, and the corresponding frequencies are 
less than 0.1 /iHz apart. Furthermore, as will be explained later on, some pulsation modes were calculated using the 
Lagrangian displacement rather than the Eulerian velocity perturbation. When comparing the two methods, differences 
on the frequencies are very small for m — and can be larger for m ^ (for example, 5uj/uj — 10"^ for the mode 
represented in Fig. 11 right panel), thereby providing yet another verification on the pulsation mode calculations. 



Overall, these tests indicate a good numerical stability both with respect to the numerical resolution and the choice 
of numerical method. The tests on the variational principle, on the other hand, show that some numerical difficulties 
remain, possibly resulting from a loss of precision on t he stellar models. F urthermore, the accuracy is not as good 
when the rotation rate approaches break-up, as shown in Reese et al. (2008b I. Before doing accurate comparisons with 
actual observations, these difficulties will need to be addressed. JNonetheless, these are not expected to change the basic 
behaviour of the pulsation modes nor the results in following sections. 



3. Uniform or nearly uniform rotation profile 

3.1. Mode classification 



As was stated above, Lignieres & Georgeot (20081 and Lignieres & Georgeot (20091 have previously shown that for 



rotating polytropic models, pulsation modes fall into the following main categories: island, chaotic, and whispering 
gallery modes. We have found that a similar classification also applies to pulsation modes in SCF models with uniform 
or mildly differential rotation (at least up to a = 0.4, which, based on Eq. (fTl), gives an equatorial rotation rate which is 
84% of the polar rotation rate). Figure [s] compares pulsation modes from both types of models. As can be seen in the 
figure, corresponding modes with an analogous geometric structure are also present in SCF models. 



3.2. Quantum numbers for island modes 



As was also the case for polytropic models, it is possible to introduce a new set of quantum numbers (n, I, m) based on 
the geometry of island modes (see lower left plot in Fig. 3]) . These quantum numbers then intervene in a new asymptotic 
formula which describes the frequency organisation of these modes: 



"Am — mfl 



fit 



(27) 



where A^, A^, and a are free parameters which depend on stellar structure. The parameter figt corres ponds to the 
rotation rate but is treated as a free parameter. This formula is quite similar to the one introduced in (Lignieres & 
Georgeot|2008 Reese et al.|2008a ) except that |m| has been replace by as this provides a slightly more accurate fit to 
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Island Chaotic Whispering gallery 

Low / — \m\ Medium I — \m\ High I — \m\ 




Figure 3. A comparison between pulsation modes in polytropic models and models based on the SCF method. The same 
three categories apply in both cases as can be seen by the analogous geometric structure. The quantum numbers h and 
I, which only apply to island modes, are the number of nodes in the directions indicated in the lower left plot. 



the frequencies. The reason why |m| had been obtained in Reese et al. (2008a I is because the formula was first derived 
for the quantum numbers (n, I, m), where a term proportional to \m\ dominates, and then adapted to (n, I, m). 

Table [T|gives the values of these parameters for selected SCF models as well as for a polytropic model. The parameters 
were calculated from a sparse frequency set and are therefore subject to some error. The ranges on the quantum numbers 
are 10 < n < 26, < ^ < 3 and —2 < m < 2. Furthermore, due to difficulties in mode identification, the parameters 
given for the most rapidly rotating models (r/ = 0.9) must be taken with caution. Nonetheless, a second calculation based 
on a more complete mode set shows that these values provide a reasonable estimate for 2 of the models (see following 
section). The true value or range of values for the rotation rate, firoai, is also provided and shows a posteriori that f2fit 
does approximately correspond to the rotation rate. 



Table 1. Parameters for the asymptotic formula Eq. (27 1 



M 

Mq 




a 


modes 


Aft 
(/.Hz) 


A,- 
Aft 


Aa 
Afi 


a 

A^ 


A^ 


^rcal 

Afi 


Afi 


poly* 


0.6 


0.0 


84 


36.7 


0.66 


0.029 


2.92 


0.827 


0.838 


0.047 


1.7 


0.7 


0.0 


40 


37.1 


0.77 


0.018 


3.52 


0.975 


0.982 


0.023 


1.8 


0.9 


0.0 


11 


33.5 


0.42 


0.011 


2.86 


1.157 


1.167 


0.038 


25.0 


0.6 


0.0 


39 


15.2 


0.79 


0.016 


3.39 


0.947 


0.969 


0.030 


25.0 


0.6 


0.2 


31 


15.1 


0.85 


0.018 


3.63 


0.944 


0.947-0.987 


0.045 


25.0 


0.6 


0.4 


31 


15.5 


0.90 


0.050 


3.37 


0.915 


0.830-0.988 


0.059 


25.0 


0.9 


0.0 


24 


12.4 


0.70 


-0.002 


3.41 


1.380 


1.387 


0.033 



Values of the different parameters from Eq. (27 1 for selected SCF models as well as for a polytropic model (first line). The first 
three columns identify the model, where r} and a come from Eq. (jl]). These parameters were based on a sparse mode set (the 
number of modes being indicated by N^adc) a nd a re therefore subject to error. The last column contains the average deviation 
between asymptotic frequencies based on Eq. ( 27 1 and the numerical frequencies. 
* Polytropic model with iV = 3, A/ = 1.7Mq and 
on the next line. 



1.84_Rq. These are also the mass and equatorial radius of the model 



As was noted in Reese et al. (|2008a l, the ratio A;~/Aij decreases for increasing rotation rate s. Using A ;/A„ = 1/2 in 



the non-rotating case 
finds a theoretical value 



l'assoul|l980|) , ai 

e of 2 for Aj-yAfi 



and the relationships between Afi, Aj and A„, A; given in Reese et al. (2008a), one 



when = 0. Since the values in Table [T| are much smaller, the "small" frequency 
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separa tion is no longer small but comparable with the large frequency separation, as was also observed in |Lignieres et aL] 
(|2006) and 'Lovekin e t al |p009| . 

In the last column, the standard deviation between the asymptotic and numerical frequencies is given. It is defined 
as follows: 



2\l/2 



modes 



i=i 



(28) 



where uj^^^"^^' are the frequencies given by the asymptotic formula and uji the numerical frequencies. Although the 
asymptotic formula captures the basic structure of the frequency spectrum (at least for low values of m), there are 
differences which are larger than observational error bars. The main causes seem to be deviations resulting from avoided 
crossings and also a slight variation of the azimuthal dependence of the frequencies with I and n (see following section) . 

3.3. High azimuthal orders 

The results presented so far were based on pulsation modes with azimuthal orders between -2 and 2. However, island 
modes also exist for high values of m as is illustrated in Fig. |4] As can be seen in the figure, high m island modes have 
an analogous structure to their low m counterparts except that they are much closer to the equator. This is similar to 
the behaviour of sectoral modes in non-rotating stars. 




M = a5.0M, 71-0.6 a- 0.0 




3>m 



7 



M=25.0Mg ?7-0.6 a -0.0 
w- 168,5/.illz m=10 



Figure 4. Two island pulsation modes, one with a low m value (left) and the other with a high m value (right). Although 
the basic structure remains the same, the mode with a high azimuthal order is concentrated much closer to the equator. 
This is analogous to what happens with sectoral modes in non-rotating stars when m increases. 



Figure [5] shows two pulsation frequency spectra with n = 15 to 20, ^ = to 1 and m = —10 to 10. The left plot is for 
a uniformly rotating model and the right one corresponds to differential rotation. The symbols represent the numerical 
frequencies and the continuous lines are a least-squares fit based on the following formula: 



n, I. m 



hAf, 



fit 



ail). 



(29) 



The term mflat has been removed fro m b oth the numerical frequencies and the fit in Fig. [5] so as to bring out their more 
subtle azimuthal dependence. In Eq. (29 1, the term with has been replaced so as to give the frequencies a hyperbolic 



dependence on m, as is visually suggested by the numerical frequencies. Besides the modification to the azimuthal term, 
the term has been removed but is compensated for by allowing the parameters Dm, fi and a to depend on I. 

Table|2]gives the values of the different parameters used to fit the frequency spectra in Fig.|5] Comparing these values 
with those in Table 1 shows reasonable agreement, provided one compares A^ and A^- with Z?„j/2/i and a{l = l)—a{l = 0), 
respectively, where the expressions Dm/2/i comes from a Taylor expansion of Eq. (29) around to = 0. Although the 
average deviations are larger in this table, Eq. 29 is a better fit to the frequencies than Eq. |27l The reason for this 



apparent contradiction is because the frequency sets used in Table |2] cover a larger range of to values. Applying Eq. |27| 
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Figures. Pulsation frequency spectrum in a uniformly (left panel) and differentially (right panel) rotating model. The 
radial order n goes from 15 (bottom) to 20 (top), and / and 1. Each hyperbola corresponds to a distinct pair {h, I). 
The symbols repr esent the numerically calculated frequencies and the continuous lines correspond to a least-squares fit 
based on Eq. (27 1. Some of the irregular features in the numerical frequencies are caused by avoided crossings. 



1/2 

to these expanded sets would yield {Sco'^) /An = 0.079 and 0.088 for the uniformly and differentially rotating models, 
respectively. 



An important difference between the two formulae, is that contrary to what is suggested by Eq. (27 1, the azimuthal 
dependence is different for / = and I — 1 modes. A likely cause is the fact that pulsation modes are closer to the 
equator at high m values. This wo uld then modify the path wh ich intervenes in the time integrals used to calculate A^, 
when working with ray dynamics (Lignieres & Georgeot 20081. As a result, a physically more relevant formula for the 



frequencies would include a A^- term which depends on m rather than an azimuthal term which depends on I. Of course, 
a quantitative calculation based on ray dynamics is needed to support this explanation. 



Table 2. Parameters for the asymptotic formula Eq. (29 1. 



Stellar parameters 



./ = 



.1 = 1 



072- 



M 



(/.Hz) 



A, 



Dfn 



a 

A^ 



2pAfi 



Drh 

Afl 



a 

A^ 



Drh 



Afl 



25 
25 



0.6 
0.6 



0.0 
0.2 



14.91 
15.01 



0.980 
0.962 



0.989 
0.953-0.993 



0.197 
0.221 



5.95 
6.04 



2.67 
2.50 



0.0166 
0.0183 



0.234 
0.244 



5.11 
4.52 



3.41 
3.52 



0.0229 
0.0270 



0.045 
0.052 



Values of the parameters from Eq. (29 1 used to fit the frequencies in Fig. [5] (i.e. the ranges on the quantum numbers are 
15 < n < 20, < Z < 1 and — 10 < m < 10). These values are similar to those in Table 111 (see text for details). The two values 
of f2rcai for the differentially rotating model (second line) are the lower and upper on the angular velocity, i.e. the equatorial 
and polar rotation rates, respectively. The parameter fiflt is twice as close to f^rcai as in Table [T] for the uniformly rotating 
model (a = 0) due to the inclusion of higher azimuthal orders. 



It is also interesting to look at what happens when h is increased to a large value. Figure [6] compares 4 sets of pulsation 
frequencies corresponding to n — 20, 40, 50 and 60. The symbols represent the numerical frequencies and the continuous 
lines a fit based on Eq. ( 29 1 . The frequencies are given in a co- rotating frame and have been shifted so that the different 
curves are at for mj ^/n — 0. Plotting the frequencies as a function of rtij \/fi rather than m causes the curves associated 
with the high order frequencies to overlap and reduces the difference between these curves and the n = 20 curve. These 
results suggest that as n goes to infinity, the azimuthal dependence of the co-rotating frequencies can be described by a 
law of the form n*/ (m/rC) where / is a function and 2^ — 5 = 1. 



3.4. An effective rotation rate 

Of particular interest is the parameter f2fit. In the uniformly rotating case, the term — mQflt represents, to first order, the 
advection of the modes by stellar rotation. The value given for fiflt in Table |2] is quite close to the true rotation rate and 
only differs by 0.93%, this difference probably resulting from the Coriolis force. In the differentially rotating case, — mfigt 
can also be interpreted as an estimate of the advection of the pulsation modes by stellar rotation. The parameter figt 
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Figure 6. Four sets of pulsation frequencies in a co-rotating 
frame in which n = 20, 40, 50 and 60. The other quantum 
numbers are I — and m — —10 to 10 for all sets. The az- 
imuthal order has been scaled by 1/ y/fi as this causes the high 
order curves to overlap and reduces the difference between these 
curves and the h = 20 curve. 



would then be an average of the rotation rate in which the weighting depends on the structure of the pulsation modes. 
We will refer to flfn as an effective rotation rate. We can then use Eq. ([T]) to calculate the position sgt where n{s) is 
equal to flat for the differentially rotating model. This is represented by the thick vertical line in Fig. [7] for the numerical 
value given in Table [2| The hashed region on either side of this line is an estimate of the error on this position using the 
difference between fl^t and fircai from the uniformly rotating model as a guide. 




M=25.0M., Tj-0.6 a 0.2 



Figure 7. The {n, I, m) = (20, 0, 10) island mode in the differ- 
entially rotating model. The vertical thick line corresponds to 
the position where fifit is equal to the local rotation rate. The 
hashed region on either side is an estimate of the error on this 
position, based on the difference between f2fit and f2reai for the 
uniformly rotating model. 



As can be seen from Fig. [7] sgt is located towards the outer regions of the star. This means that f2fit, when viewed as 
an average of the rotation profile, has a stronger weighting in these outer regions. This seems logical from the point of 
view of ray dynamics because a sound wave travelling along a ray path will spend most of its time in the outer regions 
of the star where the local sound velocity is lower. As a result, it will spend more time being advected by rotation in 
that region rather than in an inner region. Sound-travel times along ray paths have already been used to establish an 
asymptotic expression for rotational kernels of high order p-modes in spherical stars (Gough |l984 l. 

One of the best ways to confirm these ideas in a quantitative way is to apply a v ariatio nal formula which is valid for 
differential rotation. Such a formula has been established in [Lynden-Bell fc Ostriker (19671. Here we give a different and 
somewhat simpler expression which is only valid for conservative rotation profiles: 



= /" {Lo^^, + mnfpo\\ifdV + 2i [ {Lj^^, + mn)po^-{e ^i)dV 
Jv Jv 

- j Po\L?sd,{n^)dv- I / p^Nl\Q'dv+\ I llvvPiiW, 

JV JV Po^o JV JV^ 



(30) 



where ^ is the Lagrangian displacement, the displacement component perpendicular to the rotation axis and the 
displacement component in the same direction as the effective gravity. In order to apply this formula, it is necessary to 
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25M^ 77=0.6 a=0.2 



Figures. A plot of the kernel associated with the (n, I, m) — 
(20,0,10) mode shown in Fig. |7] As can be seen, the high- 
est amplitudes are reached very near the surface, just like for 
p-modes in non-rotating stars. The vertical line indicates the 
position of s'^ = f2~^ (f^eff) which can be compared with the 
position of sat as shown in Fig. [7j 



calculate pulsation modes in terms of ^ rather than v, the Eulerian velocity perturbation. This has been done, and the 
relevant pulsation equations are described in App. \K\ 

In the uniformly rotating case, the term that corresponds to the advection of pulsation modes by rotation is mf2, 
which is contained in the first integral. By analogy, we can define an effective rotation rate for the differentially rotating 
case as follows: 



Kar + mnypo\^\'dV = Kar + ™f^cff / Po^'dV (31) 

V Jv 
Solving for flag leads to the following expression: 




^oS^~^ + i:\l I {uJv.r + mnflCdV, (32) 
where 

When mH. ^ 2cj, then fioff can be approximated by Jy ^llCdV. This is similar to the first order perturbative expression 
describing the advection of modes by slow rotation, except that the kernel /C has been defined from the eigenmode in the 
rotation star. Figure [8] shows a plot of the kernel associated with the mode in Fig. [7] As can be seen in the figure, the 
highest amplitudes are reached very near the surface, just like for acoustic modes of non-rotating stars. Superimposed 
on the diagram is a vertical lines which indicates the position of s^f — fi"^ (f^cff)- This can be compared with s^t which 
is plotted in Fig. [t] As can be seen from the two figures, Sgf < Sfif This difference comes from the fact that sgt not only 
includes the advection of modes by rotation but also the effects of the Coriolis force on the mode frequencies. 

The rotation rate fioff turns out to be a very good indicator of the advection of modes by rotation. This is illustrated 

in Fig. joj which shows a comparison between I^uj-^ j — j _„i^ /2m (solid line) and (flcsim) + rieff(~'7i)) /2 (dotted 

line), for a set of modes in which the Coriolis force has been removed. The relative difference between the two is around 
10~^ — 10~^, making it difficult to distinguish the two curves. The downward trend results from the fact that the pulsation 
modes are becoming closer to the equator as \m\ increases. The dent between m = 5 and m = 6 is caused by an avoided 
crossing. Also, the curves corresponding to flcs^m) and flcg{—m) have been included. The reason why these curves are 
not identical is because modes with azimuthal orders m and —m are not identical even without the Coriolis force, because 
of the differential rotation profile. 

This natura lly leads on to the idea of applying inversion theory to probe the rotation profile using the rotational kernels 



defined in Eq. (33 1. The quantity i m ~ i -m) /2'7i is readily available from observations, once an accurate mode 

identification has been done. Furthermore, it turns out that — JyilJCdV is a very good approximation to ilcs, at 
least in the example considered above, thereby allowing the use of linear inversion theory. These kernels will, nonetheless, 
need to refined so as to include the effects of the Coriolis force. 
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Figure 9. Different measures of the effective rotation rate. The 
solid and dotted curves show an average over m a nd — to of f2off 
based on the numerical frequencies and on Eq. ( 32 1 . The two 
last curves show the separate contributions from to and — m, 
where to corresponds to retrograde modes and —to to prograde 
modes. 




cj-236.2mHz 




a- 
m = 0" 



M=25.0Me n-l.a a 1.0 
aj-267,3/xIIz m = 



Figure 10. The left figure corresponds to a chaotic mode and the right one to what appears to be a whispering gallery 
mode in SCF models with a highly differential rotation profile. No island modes are shown as they seem to have 
disappeared in the models (for low m). 



4. Highly differential rotation 

When the rotation profile becomes highly differential, the stellar structure becomes more and more deformed and the 
polar regions can, in some cases, become concave. These polar concavities result from the particular choice of rotation 
profile as expressed in Eq. ([T|) since they do not appear in models where the rotation rate increases with distance from 
the rotation axis. This deformation naturally affects the structure and organisation of pulsation modes. Figure \T0\ shows 
a chaotic and what appears to be a whispering gallery mode in models where the rotation profile is very differential. No 
island modes are shown as they seem to have disappeared. In the more distorted configurations, even whispering gallery 
modes become difficult to find. Instead, most of the modes are of a very chaotic nature. 

One way to counteract the effects of stellar distortion is to increase the azimuthal order m. Indeed, increasing the 
azimuthal order causes the pulsation modes to become closer to the equator and move away from the poles where stellar 
deformation is strongest. As can be seen in Fig. [TT] highly regular whispering gallery modes exist even in the most 
deformed configurations. Also, for models with less distortion, it is possible to find some island modes. Nonetheless, such 
modes are not likely to be visible in stars due to disk averaging effects. Therefore, if stars reach this degree of distortion, 
it will be very challenging to interpret their pulsation spectra. 



5. Conclusion 



As has been shown in this paper, results concerning pulsation modes in rapidl y rotating polytropic models can be 



gener alised to more realistic models based on the Self-Consistent Field method (Jackson et al. 2005[ MacGregor et al 
20071 provided the rotation profile is not too differential. In particular, pulsation modes fall into different categories, 
island, chaotic and whispering galler y modes, each with their own characteristic geometry, i n full agreement with previous 
calculations based on ray dynamics (Lignieres & Georgeot 2008 Lignieres & Georgeot 2009 ). The frequencies of the island 
modes obey the same type of asymptotic formula as those in polytropic models although a more careful investigation of 
their TO-dependence reveals a more complex behaviour than was previously established. This type of formula potentially 
provides a promising way of identifying pulsation modes in rapidly rotating s tars, especially at high radial orders where 
the agreement between formula and frequency is very good ( Reese et al.|200"9 l. Of course, when applying this formula to 
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Figure 11. The two figures corresponds to pulsation modes with a high azimuthal orders in models with a highly 
differential profile. The left figure corresponds to a whispering gallery mode and the right one to an island mode. As can 
be seen in these plots, pulsation modes become less chaotic with increasing azimuthal order. 

observations, one should restrict themselves to modes with low I and m values, because cancellation effects reduce the 
visibility of modes with more nodes on the surface. As a result, the approximate form given by Eq. |27| which is valid for 
low m values, is sufficiently accurate. 

A useful by-product of the asymptotic formula is an estimate of the effective rotation rate which gives the average 
advection of modes by rotation when the rotation profile is mildly differential. The obtained value indicates a stronger 
weighting near the surface, where the local sound velocity is smaller. This goes hand in hand with the intuitive picture 
based on ray dynamics that a sound wave is most advected in those regions where it spends the most time. A rigorous 
calculation based on a variational principle yields rotation kernels which confirm this picture and help provide effective 
rotation rates similar to the one obtained in the asymptotic formula, apart from the effects of the Coriolis force on the 
frequencies. These rotation kernels could then be used in inversion methods to probe the rotation profile. 

When the rotation profile is highly differential, pulsation modes tend to be predominantly chaotic, probably as a result 
of the star's geometric distortion. Increasing the azimuthal order counteracts this effect by drawing the pulsation modes 
closer to the equator thereby causing regular whispering gallery modes, and in some cases, island modes, to reappear. 
Nonetheless, such modes are not likely to be visible due to disk averaging effects thus making pulsation spectra in such 
stars difficult to interpret. 
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Appendix A: Pulsation equations based on the Lagrangian displacement 

In order to deri ve the Euler's equation in terms of the Lagrangian displacement, we begin with Eq. 13 of |Lynden-Bell Ez\ 
Ostriker (19671 and calculate its Eulerian perturbation in an inertial frame: 



Po 



Dt^ 



Pol • V {Vo ■ VVo 



(A.l) 



where ^ is the Lagrangian displacement, and other quantities have the same definitions as before. The time derivation 
operator is defined as follows: 

Dt dt 



(A.2) 



Simplifying Eq. |A.l| yields: 

[A + imnf poi + 2 [A + imn] pofl x | + po^ssds (O^) = - Vp + pgoff - PoV^-. (A.3) 

This is the same equation as w hat is used in Lovekin et al. (2009'). If one uses the following relationship between | and 
V {e.g. Christensen-Dalsgaard||2003( ) : 



dt 



(A.4) 



it is possible to show that Eq. |6]and Eq. |A.3| are equivalent. 

In terms of the coordinate system described in Section [2.3| Euler's equation takes on the following explicit form: 



= [oj + mH,]^ po 



+ 2i[uj + mn] — poC 



-poS [dsfl^) sin 6 



sine C{rg sin6l-t-rcos6') e 

^2 ^ j.2j. ^ 



Po 



p - Pod(;^ 



= [o; + mf2] Po 



n-, (rgsinO + rcos0) ^j, 
+ 2t[uj + mn] ^ ^ — — -po^ 



-PoS {dsfl'^) {rg sin( 



r cos ( 



(^siuO^ C{rgsin0 + rcos6) r 



n r , ni2 oT , ni^'^^'^''^^ cC or , f^.^Cirgsme + rcosO) g d^p 
0= [w-l-mUJ Po — E,'^ — 2i[u} + m\l\ poC,^ — 2i[u} + mill Po^ ^~ra ~ Po~ 



dep 



dgP. 



Po 



-P - Podg'i', 



nc sine 



n( {rg sin -|- r cos ( 



sin 9 



(A.5) 

(A.6) 
(A.7) 
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where ico = A. This is then supplemented by the the continuity equation, 
= p + V-M, 
the adiabatic relation, 
= p + ^-Vpo-cl {p + $-Vpo), 
and Poisson's equation, 
= A* - Ap, 

which in spheroidal geometry are: 



(sine 



(sin 9 



= (p - cIp) + ^ {d,Po - cld^Po) ^« + ^ [doPo - Cldopo) e, 



c 



y.2 _|_ ^2 



2r0 „2 



1 



2 ~^CS 



(A.8) 
(A.9) 

(A. 10) 

(A.11) 
(A.12) 
(A.13) 



Using this set of equation yields the same modes and similar frequencies as using the system of equations based on the 
Eulerian velocity perturbation v. This approach nonetheless has the advantage of yielding ^ which can be used more 
readily in the variational formula for differential rotation. 



